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We propose a new method of tangential redistribution of points by controlling 
the tangential motion in curvature driven flow of plane curves. The tangential 
velocity is designed by taking into account the shape of curves. It is based on 
^ ' combination of method of asymptotic uniform redistribution and the tangential 

velocity in crystalline curvature flow equations. 
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Introduction 



We consider an embedded and closed plane Jordan curve T which is pa- 
rameterized by a smooth function x : M/Z D [0, 1] ^ such that 
, r = Image(a;) — {x{u); u £ [0, 1]} and \dux\ > 0. We denote d^f = df/d^, 

' and \a\ — y^a.a where a.b is inner product between vectors a and b. 

^ . The unit tangent vector is T = dux/\dux\ ~ dgX, where s is the arc- 

length parameter ds = |(?ua;|du, and the unit inward normal vector is de- 
fined in such a way that det(T, AT) = 1. The curvature in the direction 
N is denoted by k. Let v be the angle of T, i.e., T = (cos sin z^) and 
N = (— sin i^, cos i^). The problem of evolution of curves is stated as fol- 
lows: For a given initial curve T'^ = Imagc(a;'') = F, find a family of curve 
{r*}t>Oj r* = {x{u,t)\ u G [0,1]} which starts from F° and evolves ac- 
cording to the normal velocity v ~ I3{x, k, v). We follow the so-called direct 
approach in which the normal velocity v is the normal component of the 
evolution equation of the position vector x: dtX = f3N-\-aT, a;(., 0) = a;"(.). 
Here a is the tangential component of the velocity vector. Note that a has 
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no effect on the sfiapc of curves, and the shape is determined by only the 
value of (3. Hence the simplest setting a = can be chosen. Dziuk^ studied 
a numerical scheme for /? = fc in this case. For general /?, however, such a 
choice of a may lead to numerical instabilities caused by undesirable con- 
centration of grid points. To obtain stable numerical computation, several 
nontrivial choices of a has been proposed. We will sketch out brief reviews 
of development of nontrivial a as follows. In Ref. 2 Kimura proposed a re- 
distribution scheme in case /? = fc by using a satisfying adu = and 
the condition (U): \dux\ = L* , where L* = Jp, ds is the total length of 
r*. Hou, Lowengrub and Shelley^ utilized (U), especially for f3 = k, and 
derived (2) with ip = 1 and ui = 0. They also commented (2) for general 
(fi with w s as a generalization of (U) without explicit statement on its 
impact on redistribution of grid points. Another scheme for /3 = fc was pro- 
posed by Deckelnick^ who used a = —du{\dux\~^). In Ref. 5 Mikula and 
the first author derived (2) with ip = 1 and w = in general frame work of 
the so-called intrinsic heat equation for (3 = (}{k,v). Moreover, in Refs. 6 
and 7, they proposed a method of asymptotically uniform redistribution for 
(} = (3{x,k,p). Besides these uniform distribution method, under the so- 
called crystalline curvature flow, grid points are distributed not uniformly 
but nicely. The second author extracted tangential velocity a = —dgP/k 
which is implicitly built in the crystalline curvature flow equation.^ The 
asymptotically uniform redistribution is quite effective and valid for a wide 
range of applications. However, from a numerical point of view, there is no 
reason to take uniform redistribution automatically. Hence a new way of 
tangential redistribution can be considered in order to take into accoimt the 
shape of curves. It is a main purpose of this paper to propose a method of 
redistribution which depends on size of curvatures. The idea is combination 
of Refs. 6 and 8. 

2. Governing equations 

According to Ref. 5, one can derive a system of PDEs governing the motion 

of curves satisfying dix = f3N + aT. Main tools are Frenet's formulae 
dsT = kN, dsN = —kT and the variable transformation 9„F = gdgf, 
where g = \dux\ is the local length. Using these tools for k = Aei{dsX, d^x) 
and (cos sin I/) = dgX, we obtain the following system of PDEs: 

dtk = d1(i + ad,k + k^(i, dtx = wd^x adsX + FN 

dtg = {-kl3 + d,a)g, dtv = d'^d^y + {a + ^1)8^^ + V^^.T (1) 
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for u G [0, 1] and t > 0, with the initial conditions and the periodic bound- 
ary condition of g, k, x for u G [0, 1] C R/Z and t) = i'{0, t) + 2n. Here 
we have assumed a particular form of the normal velocity /3 that can be 
written as: /3 = w{x, k, v)k + F{x, v). The normal velocity a will be deter- 
mined later. We have denoted by and Vx/3 partial derivatives of (3 
with respect to k^v resp., and Va;/3 = {8x^(3^0x20) for x = {xi,X2)- 

3. A curvature adjusted tangential redistribution 

The so-called relative local length r(u,<) = g{u,t)/L* plays an important 
role of controlling redistribution of grid points. Let T > be the maximum 
existence time of solution. If r{u,t) 1 holds for all u as t ^ T, then 
redistribution of grid points may become asymptotically uniform. In gen- 
eral, the uniform redistribution can stabilize numerical computation, and 
it has no effect on the shape of evolved curve. Therefore, a natural question 
arises how to redistribute dense grid points on the sub-arc where the abso- 
lute value of the curvature is large. To answer this question, we will define 
extended relative local length r^{u,t) = r{u,t)(p{k{u,t))/ {(p{k{.,t))) , where 
(F(., t)) = /pj F(s, t) ds/L* is the arc-length average of F. The function (p{k) 
is assumed to be even and positive, i.e. (p{—k) = (p{k) > and it is nonde- 
crcasing for fc > 0. In this paper, we shall use (p{k) = 1 — e + e\/l — e + efc^ 
for £ G (0,1). Note that ip{k) ^ 1 if e ^ 0+, ip{k) \k\ if e ^ 1", 
and (p{k) > (f{0) > for e G (0,1). The function ip{k) = 1 is used in 
Refs. 5-7, and the function (p{k) = \k\ is used implicitly for the crystalline 
curvature flow.® The desired asymptotic redistribution can be obtained if 
one choose a such that 9{u,t) := \nr^{u,t) holds as t ^ T. The ex- 
ponential convergence limj^r 6{u, t) = is fulfilled if one choose 6{u, t) = 



ln((e«(«^o) _ ^y-f*u,(r)dr ^ = Ki - K2di\nL* = m + n2{k0), 



where ki,K2>0 are constants. If ki = k2 = 0, then 9 is preserved. We use 
Ki > and K2 = 0, if T = oo; and use K2 > 0, if T < oo and — > as 
t ^ T. Consequently, we obtain the equation of a: 



where ip = (p{k) and ip'j^ = dk^p{k). To get unique solution a, we assume 
(a) = 0. The details can be found in the forthcoming paper. ^ 

4. Numerical scheme 

For a given A^-sided polygonal curve 7"^ = Ui!=i[^?-ii I'^t ^i'' ^^'^ polyg- 
onal curves = UiXi[^i-i) ^il (j = 1,2,...). The initial is an ap- 



ds{va) _ f if) 

if ip (if) 



+ w{t) (r-i - 1) , / = pk(3 - {dl(3 + k^(3)pl (2) 
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proximation of F" such as E for i = 1,2, . . . , N, and is an ap- 
proximation of F* at the time t = tj, where tj = J2iZo "^i time 
{to = 0) with the time step n > 0. The updated curve 'P-'+^ is deter- 
mined from by using a discretization of PDEs (1) and (2). Our scheme 
is semi-implicit scheme and the discretization is based on flowing finite vol- 
ume approach, similar to Ref. 7. Discrete quantities k^, uf and rj = \pl\, 
pI = x\ — a;f_i are constant on the ith edge »S| = [a;j_j,ccj]; and a'^ and 
xl are defined at the zth vertex x\ and are constant on the ith dual edge 
Sf* ~ [x^* , xf'^j], where xf* = {xf +xl_^)/2. The duality is defined in such 
a way that F^* = {Fj + F^_,_^)/2 for F = k,u,r, are constant on Si* , and 
a!* = (a^ + (x{_i)/2 is constant on S^. The mod(27r) periodic boundary 
condition of {vf} is — — Stt, vn+i = vi + 2'k, and other all quantities 
satisfy Fq = Fat. F,y+i = Fi. We use the notation: [drF){ = (F^^^ — F^)/r, 
(a.F)f 1 = (Ff^* - Fer)Ar\ (F^) = E;Ii Ffrf/L^ U = Y^l^rf, and 

id'snl^"'"^ = {{dsfUi - (a.F),)/rf+« where (5,F), = (f1+^ - Fflf )/r|lr ; 
and (Vi) means (i = 1, 2, . . . , A/'). 

0. {r^}, {kf}, {f^} are computed from (r, A;, z^)-' for j = and each 
resp.: {r,k,vy: 

J \ 3\ u" sgn(det(pj_i,p^+J) ^/ j\ pi 

n = Wil H = ^ ^— arccos | ^ y j , T(zy,^) = -J-. 



5^ej9 i . Prom (2) and (a) = 0, ^(fcf )ai+' = (^(fcf , V^IM*)"!^^ = 

-(Et27f *^/<P(fcf))/(Eti ri'7^(A;f )) for z = 2, 3, . . . , TV, where ^ = 
E/=2 W , and 

^ = ^ - 7^ + f - l") , = .r + .2(F/3^), 

g^ep £ {drr)l = -kiPiri+^ +c4+^ - a^^i (Vi). 

Step I {drv)\ = '^^ + (af ^* + ^ + (V./3),.r(z.f ) 

(Vi) at {xi*,ki+\iyi). 

step 5 . {drx)l = Wi{d^x)l+'+c4+\dsx)^'*+F{xliyi+^*)Niiyi+'*) (Vi), 
where Wi = wixljkf'^^* ,1^1'^^*), {d^x)l'^^ = {{dsx)t+i - {dsx)i)/rl'^^* , 
{d,x), = ixf^' - xill)/ri+\ and (d,x)i^'* = (xi^l* - xf^p/rf'*. 
Note that if w includes the term |fc| ^, it is regularized by + fc^ for 

small values of \6\ <C 1. 
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Step 6. All updated quantities {rj'^^}, {kf^^}, are re-computed by 

(r, k, 1/)-'+^. A uniform time step r is used in Step 2~5 for each j. 

5. Computational results 

We will show numerical examples. In each figures, we plot in every /xr time 
step using points, and in every 3/xr we draw polygonal curve connecting 
those points, where fi = [[T/t]/100] {[x] is integer part of x), r = 10~^, 
and T is an approximated final time which computed as follows. Let A* 
(resp. L*) be the enclosed area (resp. total length) of at time t = tj. 
In Fig. 1 and 2, T is the smallest time t satisfying A*/A'^ < 10~^, and the 
final curves are magnified by 1/V A^ . In Fig. 3, T is the smallest time t 
when max{|A7A*-^ - 1|, \L*/L*-^ - 1|} < IQ-^ is satisfied. We will give 
the number of grid points N, s G (0, 1) for <^(fc), Ki and K2 for u{t). All 
the initial points are given by uniform A/'-division of the range [0, 1] except 
Fig. 1 (right). 

Affine evolution. In case p = k^^'^, any ellipse shrinks to a point homothet- 
ically (Fig. 1). Furthermore, any embedded curve shrinks to a single point 
with the asymptotic shape is an ellipse. Fig. 2 (far left) indicates the sample 

of such evolution process. 




Fig. 1. (Left) e = and (middle) e = 0.9, of an initial ellipse with half-axes 
ratio 3 : 1, and (right) the corresponding crystalline curvature flow.''"'^ The final 
rescaled curves at f = T (bottom line). Parameters are N = 100, ki = k;2 = 100. 



Weighted curvature flow. Fig. 2 (right) indicates evolution process in case 
/? = w{v)k with the weighted function w{v) = 1 — 0.8cos(4(i/ — 7r/4)). 
Loss of convexity. In case /3 = k+F{x,v),B.n initial curve may loose its con- 
vexity in finite time for some choice of F. Fig. 3 indicates such examples 
with the test function F = —2pqsnv{q{4:x\ + x\)){—AxiS\nv + X2C0Sv), 
X = {xi,X2), p = 1.25, q = 3.0, and right two figures are in case 
F = 2pq-K cos{q'K\x\^)x.N , p = 1.956, q = 1.15. proposed by Nakamura 
et al.'^^ Note that the usual crystalline curvature flow can not realize these 
convexity-breaking phenomena. 
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Fig. 2. (Far left) and (right) are evolving curves with the initial curves used in 
Ref. 5. (Left) and (far right) are the final curves. Parameters are N = 100, e = 0.1 
and Ki = K2 = 100. 




Fig. 3. Left two figures depicts evolving curves and their final curve starting from 
the unit circle (resp. initial ellipse with half-axes ratio 2:1, right two figures). 
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